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NOMENCLATURE 


C = radial clearance of damper annulus 
D = nominal damper annulus diameter = 2R 
e = damper eccentricity 
F = X - component of damper force 

A 

Fy = Y - component of damper force 
h = damper annulus film thickness distribution 
L = damper length 

p = damper film thickness distribution 
R = nominal damper annulus radius 
t = time 

x = R0 = damper annulus circumferential coordinate 

X = X - direction radial motion coordinate 

Y = Y- direction radial motion coordinate 

z = damper annulus axial coordinate 

y = damper lubricant viscosity 

ft = frequency of vibration excitation 

SPECIAL TERMINOLOGY 

3 3 

Infinitely Long Bearing Model - axial flow is neglected 
Infinitely Short Bearing Model - ci rcumferential flow is neglected 

Driver Code - Any computer code which calls the squeeze-film damper 


force computation code 
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Section 1 
SUMMARY 

As a result of the first-year effort on this grant, a 
general purpose squeeze-film damper interactive force element 
has been developed, coded and debugged. This software package 
has been applied in nonlinear dynamic analyses of some simple 
rotor systems. 

The work completed under this first-year grant is a sig- 
nificant step in the development of strategies and add-on 
software packages which will be needed to apply available ad- 
vanced nonlinear finite-element codes (such as ADINA) to general 
engine dynamic simulation. Also, a detailed discussion is pro- 
vided of the direction of effort for the next two years. 
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Section 2 

INTRODUCTION AND BACKGROUND 
2 . 1 Engine System Dynamics 

Present day jet engine configurations have evolved pri- 
marily through a trial -and-error process involving extensive 
testing. There are many fundamental dynamic phenomena which 
take place within these engines for which basic description 
and understanding have yet to be generated. Nonetheless, they 
work well. Modern aircraft engines are typical of current 
high-technology products in which the recently acquired comput- 
ing capabilities of today are being used to better understand 
and improve what is already designed, built and operating. 

A better understanding of the basic dynamic characteristics 
of existing and new engine configurations is a prerequisite for 
producing acceptable engine efficiencies on advanced configura- 
tions (i.e. smaller rotor/stator running clearances). Also, a 
better definition of engine dynamic response would more than 
likely provide valuable information leading to reduced mainten- 
ance and overhaul costs on existing configurations. Furthermore, 
application of advanced engine dynamic simulation methods could 
potentially provide a considerable cost reduction in the develop- 
ment of new engine configurations by eliminating some of the 
tri al -and-error process done with engine hardware. 

The emergence of advanced finite element codes, such as 
NASTRAN, NONSAP, MARC and ADINA, and related algorithmic advances, 
have placed comprehensive engine system dynamic analyses within 
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reasonable reach. What remains to be done is to develop new 
component element software to properly model engine rotor/ 
stator interactive components, such as squeeze-film damper, 
within the algorithmic logic of already proven finite element 
codes. This is the major mission of this grant. 

2 . 2 The Function of Squeeze-Film Dampers 

For good reasons, aircraft engines use rolling element 
bearings exclusively. This design philosophy has, until recent 
years, deprived engines of the beneficial damping inherent in 
many other types of rotating machinery where fluid-film journal 
bearings are used. The implementation of squeeze-film dampers 
in recent engine designs has now provided engine designers with 
an effective means of vi brat ion energy dissipation. The net 
result is that the newer engines with squeeze-film dampers are 
less sensitive to residual rotor imbalance and better able to 
control vibration and transmitted force levels resulting from 
various excitation sources within the engine. 

2 . 3 Currently Available Analysis Procedures and Limitations 
The field of rotor dynamics has evolved to its present 

state primarily through the solution to problems in types of 
machinery other than aircraft engines.. In most other types of 
rotating machinery (e.g., steam turbines, centrifugal pumps and 
compressors, fans, generators, motors, etc.) the rotor can be 
adequately modelled as an Euler or Timoshenko beam^ In addition, 
the support structure holding each bearing can often be adequately 
modelled as a separate mass-damping-stiffness path to ground 
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(i.e., to the inertial frame). Also, for most purposes, bear- 
ing lubricating film dynamic properties are characterized as 
stiffness and damping elements, linearized for small vibration 
amplitudes about some static equilibrium state. It is this 
level of sophistication that has been utilized for the most 
part in rotor-dynamic analyses of aircraft engines (e.g., 

Hibner [ 2 ]). 

Present day aircraft engines are structurally far more 
complex than most other types of rotating machinery. The multi- 
shaft configuration, plus the fact that the shafts are thin 
rotating shells, not simple beams, creates unique but signifi- 
cant complicating differences between aircraft engines and other 
machinery. Also, the stator structural support at each rotor 
bearing represents anything but a separate mass-damper-stiffness 
path to an inertial frame. In fact, setting the inertial frame 
for the engine is not a simple matter when the full range of 
in-service maneuvers is realized. Dynamic paths between differ- 
ent bearings exist not only through the rotor but through several 
other paths within the non-rotating engine structure, i.e., a 
"multi-level mul ti-branch" system. As many as eight significant 
"levels" have been identified. 

The feasibility of nonlinear dynamic analyses of multi- 
bearing flexible rotors has been recently demonstrated on non- 
aircraft applications (see Adams [3]). There are highly non- 
linear dynamic effects in aircraft engines, particularly under 
large excitation forces, such as blade or disk failures, hard 
landings and foreign matter ingestion events. 
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Clearly, the field of aircraft engine dynamics is presently 
in a position where there is both a need for substantial ad- 
vances and feasible means available by which such advances can 
be accompl i shed . 

2 . 4 The Need for Time-Transient Nonlinear Dynamic Analyses 

In recent years it has become evident that an important 
class of engine dynamic phenomena can not be studied without 
accounting for the highly nonlinear forces produced at bearings, 
labyrinths and other c 1 ose- runn i ng rotor/stator clearances under 
large amplitude vibrations. In such cases, linear theory typic- 
ally predicts vibration amplitudes larger than the actual run- 
ning clearances. Furthermore, important vibratory phenomena, 
such as subharmonic resonance and motion limit cycles, are 
"filtered" out of the problem with a linear model, giving grossly 
erroneous predictions, qualitatively as well as quantitatively. 

With few exceptions, nonlinear dynamic problems must be 
solved numerically as time- trans ient responses, whether the 
sought "answer" is a steady state periodic motion or is strictly 
a transient phenomenon. The problem is mathematically categor- 
ized as an initial value problem in which the displacements and 
velocities of the complete system must all be specified at the 
beginning of the transient. From that point forward in time, 
the equations of motion are numerically integrated (known as 
"marching") as far in time as one wishes to study the system 
motions and forces. If the system is dynamically stable, the 
transient motion dies out yielding the steady state response 
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which in a system with a periodic force excitation will be a 
periodic motion. In a stable system with no time-varying force 
excitation, the transient will die out as the system comes to 
rest at one of its stable static equilibrium positions. If 
the system is unstable, the transient does not die out but con- 
tinues to grow in time unless or until some nonlinear mechanism 
in the system limits the motion to what is frequently called a 
"limit cycl e" . 

In order to study the general dynamical characteristics of 
aircraft engines, nonlinear dynamic computational schemes are 
required. The approach taken in this grant is to develop soft- 
ware packages to model engine components which are not typically 
found on dynamical structures and therefore are not already 
built into existing nonlinear finite-element structural dynamics 
computer codes. This first-year effort has concentrated on 
developing such a software package for squeeze-film bearing 
dampers . 

2 . 5 First-Year Effort, Development of Damper Element 

The main objective of the first-year effort was to develop 
a squeeze-film damper element (i.e., software package) suitable 
for implant into a general purpose nonlinear finite-element 
computer code. This objective has been met in full. Furthermore, 
workable strategies have already been developed to implant this 
damper element. Also, the damper element has been extensively 
tested on simple rotor/stator configurations under a wide variety 
of dynamic loading conditions. These results are presented in 
subsequent sections of this report. 
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Section 3 


BEARING DAMPER ELEMENT 


3 . 1 Introduction 

The bearing damper finite element code is essentially an 
interactive element to represent squeeze film dampers. That 
is, its purpose is to bridge the "gap" between structural ele- 
ments which are separated in the actual engine by a squeeze 
film damper. In its simplest version, it has an input/output 
setup as shown in Figure 1. As the bearing-damper element is 
extended to encompass more types of rotor/stator interactive 
forces (e.g., rubs, impacts, etc.) the input/output list will 
expand . 

A source listing of the bearing damper element code de- 
veloped during the first-year is given in Appendix A of this 
report. 

3 . 2 Governing Equations 

The rotor/stator interactive force generated in a bearing 
squeeze film damper is modeled using an adaptation of the class- 
ical Reynolds lubrication equation for incompressible laminar 
isoviscous films. 



( 1 ) 


z 


axial coordinate 


x 


circumferential coordinate = r0 


h = local film thickness 


-rr = instantaneous local rate of change in h 
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U = sliding velocity = Rw, typically zero in a damper 

C = radial clearance of damper annulus 

The relationship between system inertial coordinates and damper 
para meters comes through the expression for h, 3h/3x and dh/dt. 
Referring to Figure 2, these relationships are summarized as 
f ol 1 ows : 

5 = (X R - X $ )1 + (Y r - Y s ) j (2) 

• • • ^ • • A 

e = ( X R - X $ )i + ( Y R - Y $ )j (3) 

then 

h = C-e-n G = C - (X R -X $ ) cos 6 - (Y R -Y $ ) sin 6 (4) 

[< X R- X S> s1 " 0 - <VV < 5 > 

and 

^jr = - (X R -X $ ) cose- (Y r -Y $ ) sine (6) 


3-3 T ypical Configurations and Boundary Conditions (see Figures 3,4,5) 
Some engine manufacturers do not use centering springs in gen- 
eral on either military or commercial application because of fatigue. 
This can require using a tighter clearance and thus requires a 
tighter control on dimensional tolerances on annulus diameters. 

In both cases |^' << |^» i.e., axial pressure drop within annulus 
is much smaller than circumferential pressure drop. This re- 
duces the governing equation (1) to. 


d / h 3 dp , _ i o dh 

dx y dx ' dt 


(7) 


the "infinitely long" bearing equation for zero rotation (11 = 0). 


9 


Other less frequently used configurations do not employ 
end seals, in which case the "short bearing" approximation or 
its equivalent is used. In this case, the local axial end flow 
is considered to far outweigh the ci rcumferenti al flow leading 
the "short bearing" approximation. 



Actually, an improved adaptation of the short-bearing approach 
is obtained by implementing the parabolic assumption of 
O'Donoghue [9]. That is, the following approximation is made, 


which assumes an axially symmetric axial pressure distribution 
at every circumferential location. This then gives the follow- 
ing pressure field equation. 


This is actually a first-order Fourier approximation using 
the parabola as the single approximative function. 

A convergent approximation to the full two-dimensional 
Reynolds equation can be obtained, as an extension of the 
foregoing approach by O'Donoghue [9]. The number of Fourier 
terms is increased to N, resulting in N simultaneous ordinary 
differential equations. 


P ( 0 , z ) = p ( 0 , 0 ) (1 - 

L 2 


(9) 


J_ iL_ ( h 3 jLR) = i 9 — + 8 P ( Q » Q ) h 3 

M 3x 3x ; ,£ dt 0 L 2 


( 10 ) 


p(0,z) = p (0,0) cos — j- + p (0,0) cos 

1 L 2 



HD 


Substitution into the general 2-D Reynolds equation (1), expan- 
sion of the right hand side ( RHS ) by the same series, followed 
by LHR:RHS segregation by the arguments under cosine yields N 
ordinary differential equations, one for each p^(0,O). 

3 . 4 Method of Solution 

Whether the long-bearing formulation (7) or the other two 
formulations described by (10) and (11) are used, the following 
solution method is employed. It is described below as imple- 
mented for the long-bearing formulation. 

Based on 3-point central difference, 

hS 2 dh d£ ^ dh 

dx 2 dx dx M dt 

P - P 

/d£x = i + 1 i - 1 

v dxv 2Ax 


( dl£) 

dx i 


p i + i 

Ax 2 


h 3 ( >U1-^V P 1-1 , + 3h > _1 , 


dh ' P i + 1- P i-1 


Ax 


2Ax 


) = 12p 


dh. 

dt 


( 12 ) 


Rearrangi ng 


h ? 3 h ? d h . 

P f- 1 - ■ + 
i + l L . 2 + _ 

^ ; ^ / 


dh . 


. 2 h | h? 3h? dh. 


D o 


c j 


E . 
J 


r j 


Recurrsion relationship, 

CjPj + EjPj_-|+DjPj + i = Rj, form of difference equation 


Employ the form. 


P . , = A .P . + B . 
J-l J J J 


Then , 


or 


C .P . + E . (A -P . + B . ) + D .P . Al = R . 
JJ J J J J JJ + 1 J 


P.(C. + E . A . ) + E.B. + D . P . , , = R. 
J J J J y J J JJ + 1 J 


Then , 


-D, R . - E .B . 

p = i J ) p + J JL 

j ' C . + E .A . j + 1 C. + E.A. 


J J j 


J J J 


A j + 1 


B J + 1 


°j 


C.+E.A. 
J J J 

R . - E.B. 

_J 

C . + E . A . 
J J J 


(13) 


From upstream boundary condition the {A} and {B} vectors are 

determined by starting with A = 0, B = P (called forward 

2 2 1 

sweep ) . 

The downstream boundary condition is inserted at the 
beginning of the backward sweep, i.e. 


P M , - A m P m + B M 
M- 1 MM M 


Pm o - A m n P M , + B m , 
M-2 M-l M-l M-l 


P 2 A 3 P 3 + B 3 


Film rupture is handled by the following substitution. If 
P j <P vapor’ set P j = P vapor before com P utin 9 Pj. r This is 


1 2 

equivalent to the condition |-^ = 0 at the film-rupture full- 

film boundary. In the case of the 2-D convergent approach 

indicated by eqn. (11), this point- by -point test is made on 

N 

the local summation P(0.,z) = I P.(0.,z). 

J k =i K J 

The method of solution although not closed-form, is non- 
iterative. While it does entail a one-dimensional finite- 
difference scheme, it requires only a very small amount of CPU 
time and is therefore ideally suited to time transient rotor 
dynamics analyses. It has major advantages over the purely 
closed-form approximations, e.g., [10,11]. These major advantages 
are immediate account of specified-pressure boundary conditions 
at feed and drain holes of a damper. Also, the finite differ- 
ence approach easily permits account of static as well as dynam- 
ic deflections which alter the oil film gap geometry from ideal 
rigid circular shapes. 

3 . 5 Force and Force Gradients 

Forces Components on Rotor: 

e 2 

F x = - / A P cos 0dA = - LR / p ( 0 ) cos 0d0 

i 

e_ , , 

Fy = - /^p sin 0dA = - LR / 0 p(0) sin 0d0 

Stator Force Components: 

F'=-F F'=-F 

r X Y ’ r Y Y 

Force Gradients: 

3 F . 3 F . 

[C ij ] 2x2 = " ; [K ij ] 2x2 = " 3X7 

J 

AF i 3F i AF i 
ax/AXj * 3 V AX d 


(14) 


(15) 


( 16 ) 


Numerical differentiation is performed with small AX^ 
and AX. increments about instantaneous conditions. 

J 

This provides continuous updating of 

{F.}, [C..] and [K..]. See Appendix A for the computer 
J J i J 

code source listing of the completed squeeze-film element 
SQUEEZ. 


Section 4 

APPLICATION OF DAMPER ELEMENT 

4 . 1 I ntroduct i on 

For purposes of checking out the damper element code and 
to demonstrate its use, two types of computations were made 
and the results presented herein. First, a parametric study 
of damper pressure distributions was made for a variety of 
specified circular orbits, for both long-bearing and short 
bearing solutions. Second, a four-degree-of-freedom rotor- 
damper-stator model was investigated under conditions of small 
rotor unbalance through large rotor unbalance. These results 
are reported in the following sections. 

4 . 2 P ressu re Distributions for Specified Circular Orbit s 

For this series of computations the following damper annulus 
parameters were used. 

Diameter, D = 6 in. 

Length, L = 1.25 in. 

Radial clearance, C = 0.010 in. 

Lubricant viscosity, y = 1 x 10 ^ reyns 

Angle between inlet oil port and drain port, (0^ - 0 Q ) = 180° 
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Inlet oil port pressure, P . = 55 psia 

Drain port pressure, p Q = 15 psia 

Lubricant vapor pressure, p y = 1.5 psia 

Orbit angular velocity, fl = 3600 cpm (376.99 rad/sec) 

The above damper parameters are typical for modern gas turbine 
aircraft engines. A parametric study was made postulating the 

outer ring of the damper fixed and the inner ring having a 
constant-radius constant-velocity concentric orbit. Eccen- 
tricity ratios (i.e., orbit radi us/radi al clearance) from 0.05 
to 0.95 were computed, both for the long-bearing and short- 
bearing solutions (both are presently incorporated in the 
damper-element computer code). 

Circumferential center-line pressures were plotted as a 
function of circumferential position and time, for one period 
of prescribed motion. The results for the long-bearing solution 
are shown in Figure 6, and Figure 7 for the short-bearing solu- 
tion. The difference between long-bearing and short-bearing 
solutions is quite large, particularly as motion amplitudes get 
smaller. The long-bearing solution provides a considerably 
stronger damper, thus the common preference of designers to use 
end-sealed dampers. 

4 . 3 Nonlinear Dynamic Response of Simple Rotor Systems 

A simple "driver" code was written (see Appendix B for list- 
ing) which uses the damper-element code in the same manner as a 
general application with large finite-element codes. The 
"driver" code is based on a four (4) degree-of-f reedom system 


i.e., planar motion of the inner and outer damper elements. 

This then simulates a single-mass rotor connected to a single- 
mass stator via the damper element. The system analyzed is 
shown in Figure 8. The model is coded to simulate arbitrary 
rotating and/or static radial loads. Aside from demonstration 
purposes, this four (4) degree-of-freedom model has been devised 
to check against the same type of system when executed with the 
damper element implanted into the general purpose nonlinear 

finite-element code ADINA, which the University of Akron has 
purchased as its contribution to this grant. 

Note from Figure 8 that the high pressure port (i.e., 
feed port) is located on the bottom of the damper so as to 
assist "lift-off". Since centering springs are not typically 
used, they have been excluded in this example. Lift-off 
therefore requires some amount of vibration to overcome the 
dead weight load. Rotating unbalance loads of 100, 200, 300, 

500 and 1000 lbs were run with = 150 rad/sec. Orbital plots 
were made showing rotor and stator total motion on one plot and 
rotor-rel ati ve-to-stator motion on a second plot. The plotted 
results are shown in Figure 9 through 13. 

For a 100 lb rotating load (Figure 9) the motions shown are 
for a 20 cycle transient from time = 0. The rotor and stator 
each show close to the same motion, and their relative motion 
is small, with the rotor barely "lifting off". The relative 
orbit is essentially oscillatory. However, when the rotating 
load is increased to 200 lbs, (Figure 10), the relative orbital 


motion shows the beginnings of orbital motion, i.e., a "cresent 
moon" shape as measured by numerous investigators. Further 
increase in magnitude of the rotating load to 300 lbs (Figure 11) 
shows a well defined steady-state total motion as well as rela- 
tive motion. Note that with a 300 lbs rotating load, the rela- 
tive ( rotor-to-stator) orbit is still small in comparison to the 
radial damper clearance and confined to the region of the bottom 
of the damper. However, an increase of rotating load magnitude 
to 500 lbs causes a considerable change to the relative orbit 

(Figure 12). Notice now that the relative motion of the rotor 
with respect to the stator fills a major portion of the clear- 
ance circle. Further increase of rotating load magnitude to 
1000 lbs (Figure 13) simply causes the steady-state relative 
orbit to expand and fill even more of the damper clearance 
circle. 

Section 5 

FINITE ELEMENT IMPLANT STRATEGY 

The previous sections gave a thorough discussion of the 
development of the interactive squeeze film bearing element. 

This section will outline ongoing efforts aimed at incorporat- 
ing these elements into the finite element procedure. In this 
context, the discussion will be organized into several main 
parts , namely : 

i) Choice of FE code used for initial implantation; 

ii) Overall solution strategy; and, 

iii) Solution algorithms employed. 
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5 . 1 C hoice of FE Code Used for Initial Implantation 

Before discussing the choice of FE code adopted, it is 
worthwhile to briefly overview various of the salient features 
associated with rotor-beari ng-stator modelling. To organize 
our thoughts, we consider them in two main phases, namely: 

i) Normal operating conditions; and, 
ii) Abnormal operating conditions. 

For normal situations, since the clearance between the blade 
and shroud and the various engine seals are quite small, the 
overall kinematic description can be cha racteri zed by small 
strains superposed on an initially small field [4]. Because 
of this, except for local zones, the overall structural mater- 
ial characterization can be considered essentially Hookean in 

nature. In this context, the structural modelling of the engine 
can be considered essentially linear in nature. Regardless of 
this though, as has been seen from the discussion in the pre- 
vious sections, even small unbalance loads can initiate highly 
nonlinear interactive forces in the squeeze film bearings. 
Because of such nonlinearity, under normal operating conditions 
the rotor- beari ng-stator system can be modelled as a partitioned 
system wherein the structural components are linear while the 
bearings are nonlinear. 

For abnormal operating conditions, the rotor excursions 
are on the order of the various blade and seal clearances. In 
this context, due to the relative smallness of such clearances, 
the deformation process can be characterized by at most small 


) 
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strain moderate rotations superposed on small initial fields 

[12]. Since such kinematic excursions are still deemed small, 

★ 

except for local events , the global structural material be- 
havior can still be considered Hookean. Because of this, the 
structural modelling of the engine can be assumed kinematically 
nonlinear. Obviously, during a fatal event both kinematic and 
massive material nonlinearity are evidenced during structural 
col lapse. 

In the context of the foregoing, it is of utmost importance 
that the FE test code chosen have adequate nonlinear element 
substructural capabilities to allow for the proper partitioning 

into linear and nonlinear element groups. This obviously 
enables more efficient running characteristics. Together with 
the partitioning capabilities, the code should also have an 
efficient updating architecture. As this feature is typically 
the heart of any nonlinear solution strategy, it is an abso- 
lutely essential characteristic. In addition to the foregoing 
features, the code chosen to test the bearing element should have 
[5]: 

i) Accessible program architecture; 

ii) Efficient running characteristics; and, 

iii) Flexible algorithmic options. 


* blade impacts, creep/fracture of blades 


Since many general purpose codes such as NASTRAN, STRUDL, 
FESAP, etc. are essentially linear with grafted nonlinear capa- 
bilities, they tend to have a less efficient/flexible program 
architecture. Because of this, our attention must turn to codes 
such as ADINA, ANSYS, MARC, etc. Since ANSYS and MARC have 
somewhat inaccessible program archi tecture , the ADINA program 
was chosen to check out the "bearing element implant". This 
follows since ADINA has the requisite combinations of capabil- 
ities, namely [6]. 

i) Nonlinear element partitioning feature; 

ii) Efficient updating architecture; 

i i i ) Flexible algorithmic options; 

iv) Accessible program architecture; and, 

v) Efficient running characteristics. 

5 . 2 Overall Solution Strategy 

The initial approach taken has been to implant the bearing 
element directly into the ADINA architecture so that direct 
numerical time integration algorithms can be employed to generate 
the transient rotor stator solution. To simplify the discussion, 
the presentation will be organized into several main areas, 
namely : 

i) Element architecture; 

ii) Overall FE code architecture; and, 

iii) Solution methology 
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The overall architecture of the bearing element is being 
structured to have several main options, namely: 

i ) Initial I / 0 ; 

ii) Interactive I/O; 

iii) Generalized stiffness and damping connectivity; and, 

iv) Generalized element library. 

The initial I/O options involve a one-time input of various pre- 
selected parameters including such categories as: 

i) Geometric configuration; 

ii) Material properties; 

iii) Element selection; and, 

iv) Required element connectivities. 

Each of these categories are in turn broken down into several 
different items, for instance: 

1. Geometric Configuration 

i) Inner and outer damper radii 

i i ) Bearing length 

iii) Orientation of oil feed grooves 

iv) Structural clearances 

v) Placement of roller bearings 

2. Material Properties 

i ) Oil properti es 

ii) Temperature dependence 

iii) Roller bearing force deflection characteristics 
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3. Element Selection 
i) Short bearing 

i i ) Infinite bearing 
iii) Roller bearing characterization 
iv) Rub/ impact 

4. Required Element Connectivities 

Several of the foregoing parameters are being coded to be 
interactively redefined depending on the nature and level of 
excitation for example, temperature and structural clearances 
fall into this category. Additionally, such field variables as 
film forces as well as the instantaneous tangent stiffness and 
damping matrices are being coded so as to be interactively re- 
defined. Such parameters are up dated depending on the nature 
of the interactively calculated position and velocity histories. 
In this context, the various interactive field quantities now 
being coded into the bearing implant associated with the ADINA 
code consist of: 

1. Velocity differential developed across the squeeze film; 

2. Positional differential developed across the squeeze film; 

3. Interactive force field developed; 

4. Tangent stiffness matrix developed by suqeeze film; 

5. Tangent damping matrix developed by squeeze film; and, 

6. Tangent stiffness of roller bearing. 

To generalize the capability of the "bearing implant", 
the initial and interactive I/O modes of data transfer are 
being developed so as to admit fairly extensive structural con- 
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figurations. This includes the possibility of accessing the 
entire array of structural elements inherent to ADINA together 
with the various constitutive models including 

1 . Hoo kean 

2 . Plasticity 

3. Temperature dependent properties 

4. Mooney Rivilin [3 ], etc. 

The overall architecture of the bearing implant is defined in 
Figures 14 and 15. As can be seen in Figure 14, the bearing 
implant is being imbedded in a buffer routine which will serve 
primarily as a link between the various data transfer modes of 
ADINA, namely: 

i) Common blocking (dynamic form); 
ii) Subroutine parameter lists; 
i i i ) Disk I / 0 . 

The buffer routine will also serve to convert the interactive 
information into the appropriate partitioned form for assembly 
into the mainstream of data flow inherent to ADINA. Namely, 

the tangent stiffness and damping matrices together with the 
interactive forces will be assembled into the proper locations 
in their global counterparts. This is currently being programmed 
into both the in core and out of core storage mode options in- 
herent to ADINA. 

Additionally, the buffer routine will be programmed to con- 
tain a degree of adaptive updating which will enable a more 
accurate calculation of the tangent stiffness and damping matrices. 
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Specifically, since the stiffness and damping matrices are calcu- 
lated by admitting a perturbation in the position and velocity 
fields of a given state, care must be taken to insure that the 
perturbation is neither too small nor too large. In the case 
that the perturbation is too large then the stiffness calculated 
will act more like a secant stiffness and hence be inaccurate. 

If too small, then roundoff error may be introduced into the 
calculations. To circumvent this difficulty, the current and 
past fields are compared. If the percentage changes are deemed 
too large/small, then the levels of perturbation introduced 
can be either contracted or expanded to insure proper evalua- 
tion of the tangent matrices. 

While the structure of the buffer will be somewhat depend- 
ent on the ADINA architecture [7] the main core of the beam implant 
will be more or less code independent. The actual flow of data 
into the core of the implant is achieved by subroutine argument 
lists. Figure 15 defines the overall flow of control within 
the core program of the bearing implant. The architecture of 
the core program is being made flexible enough to admit new 

options as they become available. 

Based on the foregoing bearing element implant, the 
architecture of the overall FE code is defined in Figure 16. 

As can be seen, the overall flow of control is broken into 
several major steps, namely: 

1. Initial I/O, including: 

i) Structural information 
ii) Bearing information 
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iii) Boundary conditions 
iv) Applied load 

v) Element connectivities 

2. Structural element generation, including: 
i ) Linear el ements 

ii) Partitioned assembly of linear elements 
iii) Nonlinear stiffness update loop with partitioned 
assembly of nonlinear structural elements 

3. Bearing element generation, including: 

i) Tangent stiffness and damping matrix generation 

ii) Development of right-hand side loads 

iii) Partitioned assembly 

4. External load generation 

5. Integration algorithm, including: 

i) "Stiffness" inversion 

ii) Implicit integration 

a) Newmark 

b ) Wilson 

iii) Explicit integration 

a) Central difference 

6. Convergence checks 

i) Norm test of out of balance loads and nodal displace- 
ments 

ii) Higher order checks 


7. Clearance checks 
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8. Adaptive Strategies, including [7]: 


i) 

Structural stiffness updating 

ii) 

Bearing 

stiffness updating 

1 1 1 ) 

Choice 

of perturbation size 

i v) 

Choice 

of integration algorithm 

v) 

Choice 

of time step size 

vi) 

Preferential partitioned updating, 


A simplified view of the actual flow of control is given 
in Figure 17. This figure includes both the linear and nonlinear 
structural loops. Currently such modifications are being in- 
serted into the ADINA architecture. 

5 . 3 Solution Algorithms 

As noted earlier, having developed the "bearing element", 
the current thrust is to implant the element into ADINA wherein 
direct numerical integration will be employed to generate the 
transient solution. In this context, several types of integra- 
tion operators are being incorporated into the coding. In 
particular various versions of the following operators are 

being considered: 

1 . Newmark [9] 

2. Wilson [10] 

3. Houbolt [11] 

4. Central difference [11] 

5. Hughes [12] 

6. Felippa, Park, etc. [13] 
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Section 6 

DISCUSSION - DIRECTIONS OF FUTURE WORK 

In view of the modelling deficiencies noted earlier, a 
more direct way of handling the structural aspects of the rotor- 
beari ng-support ( RBS ) system is necessary if a proper transient/ 
steady state model is to be developed for jet engines. In this 
direction, it appears that the finite element (FE) method is the 
requisite modelling approach for such problems. This follows 
from the fact that its inherent capabilities include the follow- 
ing features: 

i) The FE procedure has the capability to handle mul ti -branch/ 
level structure in a more direct and efficient manner than 
flexibility approaches; 

ii) The approach is well suited to handle nonlinearities due to 

a) kinematic and kinetics associated with the structure 
[14]; 

b) various types of boundary and constraint conditions 
[14], and ; 

c) material characterization [14,15]. 

iii) A body of established and proven algorithms which can handl 
various types of nonlinearities has evolved; this includes 
both the capability to handle static [14,15] as well as 
transient situations [14,16]; 

iv) Modelling of overall RBS systems more direct as extensive 
element libraries are currently available; this includes 
beam, plate, shell, 2-D, as well as 3-D elements [15]; 
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v) Algorithmic adaptability. 

Currently available general purpose codes such as NASTRAN, 
MARC, ANSYS, ARGUS, ADINA, ASKAII, NEPSAP, FESAP, SAPVI all have 
most of the foregoing items implemented as user features [17]. 
Although these codes possess the required degree of generality 
to model the structural aspects of jet engine rotor-stator 
structure, what is currently lacking are interactive "bearing 
type elements" and the overall algorithmic strategies to handle 
conservative/nonconservative interactive type forces. In turbine 
engine, such fields are generated in the squeeze film dampers 
and labyrinth seals and during rub-impact events. 

In addition to the foregoing modelling difficulties, there 
is also a need to better quantify the effects of such factors as: 

i) Rotor/stator static de-centering forces generated via: 

a) manufacturing tolerances 

b) thermal warps 

c) high "g" forces 

d) in service damage and wear, 

ii) Degree of structural nonlinearity encountered, 

iii) In service dynamic phenomena (rubs, impacts, etc.). 

6 . 1 Compatibility With Proven Finite Element Codes 

As noted earlier, while currently available FE codes possess 
the requisite generality to handle the structural aspects of RBS 
system modelling, no provisions are currently available to model 
the conservative/nonconservative effects of squeeze film damp- 
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ers, seals, rubs, impacts, etc. In view of this, future efforts 
will be given to developing a variety of special purpose 
"bearing elements" which can model such rotor/stator interactive 
force fields. These "elements" will be developed so as to be 
both algorithmically as well as architecturally compatible with 
proven FE codes. In this direction, it appears that codes such 
as ADINA would be the most likely software condidates about 
which such a development should be configured. This follows from 
the twofold fact that such codes have the following: 

i) Extensive and well proven dynamic/element capacity, 

ii) An architecture developed to allow the user to modify the 
overall algorithmic flow of a given solution loop. 

6 . 2 Preliminary Engine Dynamics Analyses 

The computational schemes ultimately implemented to track 
engine dynamic response will have to function properly over a 
wide spectrum of motion frequency and a wide range of nonlinear- 
ities. The development of computationally reliable interactive 
elements, such as the bearing/damper element, will therefore 
require a simplified engine dynamics analysis, using available 
rotor- dynamics computer codes, to realistically assess potential 
computational difficulties. For example, specifying the outer 
envelope or limits of the bea r i ng/ damper element must be pre- 
dicated on a correct understanding of relative rigidities and 
dynamic participation of individual components in and around the 
bearing. These analyses include the following: (i) linear un- 

balance forced response, (ii) linear nonsynchronous forced 
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response, (iii) linear self-excited instability analysis, and 
(iv) simplified ti me- tra ns i ent nonlinear analysis. 

6 . 3 Interactive Elements for Labyrinth Seals 

The typical jet engine configuration contains several laby- 
rinth seals. The flow field within these seals results from the 
combined effects of rotation and pressure-gradient induced axial 
through flow. Depending upon the design parameters of a laby- 
rinth seal, either a centering or decentering static radial 
force can be produced on the rotor. Likewise, the mechanical 
impedence (stiffness, damping and virtual mass) between rotor 
and stator at the labyrinth seal is a strong function of design 
details. Carefully conducted experiments by Wright [18] have 
recently shown that the labyrinth aerodynamic forces can be 
either stabilizing (positive damping) or destabilizing (negative 
damping) depending upon the direction of entering flow pre-swirl 
and the direction of rotor whirl. 

The full importance of labyrinth seals to total engine 
dynamic analyses is therefore not confined only to the potential 
for rotor/stator rubs and impacts under high vibration levels. 

A realistic simulation of engine dynamic phenomena, linear as 
well as nonlinear, must therefore include a comprehensive mathe- 
matical model for the labyrinth seals which are located through- 
out the engine. The development and implementation of a laby- 
rinth-seal interactive element is therefore important future work. 
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6 . 4 Rotor-Stator Rub/Impact Elements 

While significant efforts have been given to developing 
codes which can handle the impact behavior of compressor blades, 
no work is currently available on modelling rotor-stator rub- 
impact events. Due to the structural flexibility and close tol- 
erances inherent to gas turbine engines, such phenomena must 
undoubtably play an important role in defining the transient/ 
steady state behavior during moderate and large excursion situa- 
tions. Because of this, in addition to developing "bearing 
elements" some attention must be given to FE modelling the rotor- 
stator rub- impact events occurring in the labyrinth seals, and 
blade-case zone. Such "rub-impact elements" will have to be 
capabl e of : 

i) Tracking the appropriate rotor-stator clearances 

ii) Model impact-detachment mechanisms 

iii) Model traction and kinematic constraints generated during 
rubbing 

iv) Properly model energy losses occurring during such events 

6 . 5 Rotor/Stator Static Radial Offsets Loads 

The stiffness and damping characteristics of fluid film 
bearings are highly dependent upon their static centering or 
de-centering loads. Clearly, the squeeze-film dynamic forces 
will change considerably as static radial load is applied at 
the bearing since a shift of equilibrium eccentricity position 
will occur. Prominent sources of static radial offset loads 
result from each of the following: 
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i) Manufacturing and assembly tolerances 
ii) Thermal distortions 
i ii ) High g- force 
iv) In-service damage and wear 
v) Aerodynamic forces 

An evaluation of these static radial offset loads is in 
itself a major effort. However, a realistic computer-simulation 
of various engine dynamic phenomena can not be accomplished 
without a successful effort to determine the static inter- 
active forces between rotors and stator. 

6 . 6 Structural Non 1 i nea r i t i es 

In addition to modelling nonlinearities induced by the 
rotor/stator interactive force fields, purely structural effects 
may also be encountered. Such structural nonlinearities fall 
into two main categories: 

i) Kinematic and kinetic (geometric) [14] 
ii) Material characterization; plasticity, viscoplasticity 

The ki nema t i c- ki net i c characterization itself falls into 
three main categories, namely small deflections, small stra ins- 
large rotations and moderate/large strains. Apart from highly 
localized events such as impact-rub zones, the most prevalent 
geometric modes will most probably be typified by small de- 
flection or at most small strain-moderate rotation character- 
izations. 

For localized rotor/stator rub-impact zones, in addition 
to interactive traction fields and surface machining, potential 
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plastic flow and moderate straining may occur. Beyond inducing 
outright failure, such localized fields may have a significant 
enough effect on the geometric configuration as to cause engine 
imbalance. 

Because of the foregoing, the potential existence of kine- 
matic-kinetic and constitutive nonlinearity must be accommodated 
in the overall model. Since the effects of such nonlinearity 
are fairly well localized, a partitioned linear/nonlinear approach 
should be employed for the finite element model. 

6 . 7 Dynamic Loads 

Emergency modes of operation, such as occur with blade 
failure, hard landings and foreign matter ingestion events, 
will require a comprehensive investigation to identify and model 
the resulting dynamic input loads to the engine system. Some 
worthwhile information could be obtained from a comprehensive 
engine dynamics simulation, even with postulated high amplitude 
dynamic input loads, such as the relative endurance of two dif- 
ferent engine configurations. However, real-event simulation 
will require an accurate prior appraisal of dynamic load inputs 
to the engine system which result from identifiable emergency 
operating modes. The effort required to determine reliable 
estimates of dynamic input force time or frequency signature 
could be substantial. 

As noted earlier, the engine structure must survive a 
rather severe operating environment. In addition to extreme 
thermal and aerodynamic loads, the RBS system may be subject to: 
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i) Transient and steady state inbalance loads 

ii) Rotor-stator rub-impacts 

iii) Rotor-stator decentering forces 

Such events are generally caused by a combination of the follow- 
ing broad categories of factors: 

i ) Blade erosion 

ii) Blade-disk-seal failure 

iii) Thermal warps of rotor-stator structure due to ratchet- 
ing and creep 

iv) Misalignments due to manufacturing tolerances 

v) High "g" loads due to maneuvering 

In view of the foregoing, future analytical modelling of 
RBS systems inherent to engines must employ proven computa- 
tional schemes which possess the capability to handle as wide 
a cross-section of the loading environment as possible. As the 
time history of such loading events covers a wide range of time 
scale, the overall approach must also possess a high degree of 
algorithmic adaptability so as to accommodate both explicit and 
implicit integration schemes [16]. This is of potential import- 
ance since such schemes have been found to have varying degrees 
of success over various time scales [16]. 

6 . 8 Simulation of In-Service Dynamic Phenomena 

The direction of future work outlined here will represent 
a major advancement in the state-of-the-art of engine system 
dynamic analysis. Proper account of structural complexities, 
various rotor/stator interactive forces (static and dynamic), 
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important nonlinearities, aerodynamic forces and well defined 
dynamic load inputs will provide a greatly expanded scope in 
the types of engine dynamics phenomena that could be studied. 
Engine configuration improvement studies which are impractical 
to accomplish through testing can then be pursued through 
systematic dynamic simulation studies. 

Also, a better understanding of the dynamic behavior of 
existing engine configurations can provide valuable information 
leading to major reduction in engine maintenance and overhaul 
costs. Engine dynamic behavior is becoming progressively more 
important as efficiency improvement considerations push rotor/ 
stator running clearances progressively smaller. A realistic 
evaluation of potential engine reliability degredation result- 
ing from smaller rotor/stator running clearances demands the 
high level of dynamic system simulation described here. 

Section 7 
CONCLUSIONS 

General engine dynamic analyses which properly account for 
rotor- to-stator and rotor-to-rotor interactive forces can be 
approached through the use of available general purpose nonlinear 
finite-element computer codes. Interactive forces originating 
at bearing squeeze-film dampers and rub- impact events are, how- 
ever, not available with general purpose codes at this time. The 
work described herein shows the viability of using general purpose 
f ini te-el ement codes for engine dynamic analysis. Also, the four- 
degree-of-f reedom example model demonstrates the use of the 
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squeeze-film damper code developed in this work. Results with 
this demonstration model are consistent with the results of other 
investigators of nonlinear squeeze-film damper dynamics. 
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Damper Element Fortran Listing 
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LI ST UTILITY 


C 

C 


SUBROUTINE S CUEEZ ( AC ♦ AL , AC . A ISC , ATH 1 , ATH 2 , AE 1 .AB2,U,V, 

1 UD T . VDT »UB* VBi liBTi VBT ( AK1 1 » A K 1 2. AK22. AC 1 1 . AC 22 , F 1 , F2 . NGRIA, NSOLA , 

2 NPORA ,KAFK. KAFC.NFILA, PVAA) 

NONLINEAR TIME-TRANSIENT SQUEEZE-FILM DAMPER INTERACTIVE ELEMENT 
THIS CODE COMPUTES INSTANTANEOUS FORCE VECTOR AND ITS 
SPATIAL GRADIENTS. I.E., THE TANGENT STIFFNESS AND DAMPING 
MATRICES. 


NOMENCLATURE 


INPUT 


BD=NCMINAL DAMPER ANNULUS DIAMETER! IN) 

3L=NCMINAL D AMPER ANNULUS LENGTH(IN) 

B C=0 AMPER ANNULUS RADIAL CLE A RANCE ( I N ) 

VISC=DAMPER LUBRICANT V I SC OS I TY ( REYNS ) 

P V AP =F I L M RUPTURE P R E S SURE ( P S I A ) 

THT( 1 )=DOS IT ION ANGLE OF LUBRICANT PORT-1 (DEG) 

THT ( 2 )-= D OSIT ION ANGLE OF LUBRICANT PORT-2! DEG) 

PE ! 1 ) = SPECIF 1 ED BOUNDARY PRESSURE AT PORT-l!PSIA) 

PE! 2 ) = SPEC IF IED BOUNDARY PRESSURE AT PORT — 2 ! PS I A ) 

N G R I D= NUMB ER OF FINITE-DIFFERENCE GRID POINTS PER ’ DAMPER AR C ! ODD) 
-1 

NSGLN= l .LONG-BEAR I NG SOLUTION USED 

-2. SHORT-BEARING <PARA30L IC ) SOLUTION USED 
= 3 * FOUR IER— SERIES 2-D CCNVERGENT SOLUTION USED 
NPORT=NUMEER OF LUBRICANT PORTSlO.l OR 2) 

IF NPORT-O, JO I NED-BOUNDARY CONDITION IS USED 
NFILM=NUMBER OF IDENTICAL ANNULI FOR THE DAMPER 
KUFK=0 . ST IFFNESS MATRIX NOT COMPUTED 


AnnnAnnnonnnr>nnnnnnAnArtnAn nnnAnnnonnnnnnnnn^nonn 
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LIST U^ I LI TV 

*UFK=1 ,ST IFFNESS MATRIX COMPLTED 
K.0FC = 0 .DAMPING MATRIX NOT COMPUTED 
KOFC=l .DAMPING MATRIX COMPUTED 

X = X- INERT I AL COORDINATE OF DAMPER INSIDE SURFACE CENTER-L INE ( I N ) 
Y=Y— INERT I A L COORDINATE OF DAMPER INSIDE SURFACE CENTER-L INE. { 1 N ) 
XDT= X- IN ERT I AL VELOCITY OF INSIDE SURFACE CENTER-L I NE ( IN/SEC ) 
VOT=Y- INERT1 AL VELOCITY OF INSIDE SURFACE CENTER-L I NE ( I N/SEC ) 

XB=X- INERT I AL CUORDINATE OF CAMPER OUTSIDE SURFACE CENTER-L I NE ( I N ) 
YB=Y— INERT IAL COORDINATE OF CAMPER OUTSIDE SURFACE CENT ER-L 1 NE < 1 N ) 
XET=X- INERTI AL VELOCITY OF DAMPER OUTSIDE SURFACE CENTER-L I NE< IN ) 
YBT= Y- INERTI AL VELOCITY OF DAMPER OUTSIDE SURFACE CENTCR-L I NE ( I N ) 

/ 

OUTPUT 

Fl-X-FORCE COMPONENT CN INSICE DAMPER SURFACE(LBS) 

F 2=Y— F OR CE COMPONENT CN INSICE DAMPER SURFACE ( LBS ) 

FORCE COMPONENTS ON OUTSIDE CAMPER SURFACE ARE EQUAL OUT OPPOSITE 
FI AND F2 BECAUSE FLUID INERTIA EFFECTS ARE NEGLECTED 
SYMMETRIC PORTION OF STIFFNESS MATRIX: 

AK1 1=KXX(L3S/IN) 

AK1 2=KXY( LBS/ IN ) = KYX 
AK22=KYY ( LBS /IN) 

DIAGONAL PORTION OF DAMPING MATRIX 
AC 1 1 = CXX(Lct*5EC/ IN) 

AC 2 2 = C Y Y ( Ld* SEC/ IN) 

IMPLICIT REAL * d (A-F.O-Z) 

CC MM ON / INPUT 1/dD.dL .BC.VISC .THT (2) »-PE ( 2 ) . P V AP 
CCMMCN/I NPUT2/NGR1D .NSOLN. NFCRT ,NF ILM 
CCMML'N/ COORD /X ♦ Y , XDT » YDT ,Xd . Y3.XBT . YBT 

C > 

COMMON /F ILM/ Th( 1 0 1 ) , h ( 1 0 1 ) . QhDX ( 1 01 ) • OHO T (101 ) , 3TH( 10 1 ) . CTH( 10 1). 


non nr\ n n n nn 
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LIST UTILITY 
C 


C 

C 

C 

C 

C 

C 

C 


1DXD( 2 ) * ALF A ( 2 ) 

CCMMON/WQRK/FNGDM 1 * KB * KGUN T 
CCMMON/I NC/HMIN »V EL * C EL S ♦ D EL S T 

DIMENSION A ( l 01 i »B( 1 01 ) 4 C( 10 1 ) *E( 101 ) * RH < 101) , P ( 1 0 1 ) « ARG 1(10 1) , 
1 ARG2 ( 10 1 ) ,D < 101 ) 


ALLOCATE INPU1 

BO 

- 

AD 

6 L 

= 

AL 

3C 

= 

AC 

V ISC 


A ISC 

T h T ( 1 ) 

= 

ATH 1 

THT< 2 ) 

•— 

A TF2 

PE ( 1 ) 


A B 1 

p e ( 2 ) 

= 

A 02 

P V AP 


PVAA 

X 


U 

Y 


V 

XDT 

= 

UDT 

YDT 

sr 

VDT 

xe 

•= 

LB 

YE 


VB 

XET 


UBT 

YET 

= 

VBT 

NGRID 


NGR I A 

NSGLN 


NSOLA 

NF OR T 

= 

NFOR A 

KCFK 

.= 

KAFK 

K OF C 

= 

K AFC 

NF ILM 

= 

NF IL A 

P I 


2 • 1 A 1 5 

*R ITE 

( 5 

*2500 ) 


3 

1 

2 


' ' * ’ ~ T » ■ 1 ' I ^ I X * * r \ — J i n I ur\ ♦ f I — OIMI 

5 X » • FOTOR XCT • , 

5X, ‘ROTOR YDT • , 4X, • STATOR XD T • ,5X , • ST AT OR YDT'.//, 
5X,8(2X,D13.6)J 


*k 1TE ( e , 1 S9C ) 

1999 FORMAT 


4X» 'BEARING ELEMENT INFORMATION' , / / , 4 X » 

BL BC VISC ' 

T H 1 TH2 PBI PB2 

NS NP NF KK KC ' ) 

WRITE(6, 10)bD,8L,BC,VISC,PVAF,THT(l) , T HT < 2 ) ,PB { 1 ) , PB ( 2 ) , 



NG 


10 FORMAT 
SET UP 
KCUNT 


NGRID.NSOLN.NPORT *NF I LM . KOFK ,KOFC 
(9C1 1 .4,415.214 ) 


= I 


aaa aaa a a 
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LI UTILITY 


C 

c 

c 

c 

c 

c 


c 

c 

c 


MTEST = NGRID + 1 

NTEST - MTEST/2 

KTEST = 2* NTEST 

I F { KTEST ,NE. MTEST ) NG R I D=NGfi I O- 1 

FNGDM1 = FLOAT ( NGR I D** 1 ) 

THT(l) = THT( 1 ) *P I / 1 80 » 


IFCNPGRT.LT.2) GO TO 20 


T HT { 2 ) = THT(2)*PI/ l£C. 


IF ( TFT ( 2 ) .LT • THT ( 1) ) THTI 2 ) = THT (2 ) +2 .*PI 


ALF A( 1 ) 
ALFA I 2 > 
DXD ( 1 ) 

O X D ( 2 ) 


THT( 2 )-THT( 1 ) 
2.+PI-ALF A ( 1) 
O.S*BD*ALFA(l )/FNGCMl 
0 • 5*BD* ALF A ( 2 ) /F NG CM 1 


GO TO 40 


20 DXD ( 1 ) 
ALF A < 1 ) 


ED * P I /FNGDM 1 
2.*PI 


40 CONTINUE 


WRITE (6, 10)0X0(1 ) » D XD ( 2 ) 

A K XX = 0.0 

AK.XY = 0.0 

A KY X = 0.0 

AKYY = 0.0 
A C XX = 0.0 
AC XY = 0.0 
ACYX = 0.0 

ACYY = 0.0 
FX =0.0 
FY = 0.0 

ERANCH ACCORDING TO SOLUTION DESIGNATED 

60 CONTINUE 

CALL INCRNT 

GO TG ( 100. 1 C0.300) .NSOLN 
100 CONTINUE 

SodvE FOR SQUEEZE FILM PRESSURE DISTRIBUTION 


DO 190 KB= 1 . NPCRT 
P( 1 ) =PB( K3) 

JFCK6.EQ.1) PC NGRID )=P8( 2) 
IFCKB.EQ.2) P( NGR I D )=PD ( 1 ) 
D X = DXD ( KE ) 

CALL DFILM 
A ( 2 ) = 0 .0 
B { 2 ) = PCI) 

DC 110 K=2» NGRID 
C0EF1 = HC K) ** 3/DX* *2 


non 


43 


LI $T 


1 1 0 

120 

1 30 
C 

1 35 
1 40 


150 

160 


190 

500 

520 


UTILITY 


COEF2 = ( 1 . 5*H( K ) **2 ) *DHDX ( K )/D X 

C ( K ) =-2.*COEFl 

0 ( K ) = COEF 1 +COEF2 

E(K) = CCEF1-COEF2 

RH(K ) = 1 2 .* V I 5C +OHDT ( K ) 

I F ( NSOLN .EC . 1 ) GO TO 110 

C(K) = C ( K ) - ( 8 . *H( K ) ** 3 )/BL * *2 

CONTINUE 

NGRD = NGQID-1 

DC 120 K = 2 * NGRD 

FOCTR - C ( K ) +5 ( K ) * A ( K ) 

A( K+ 1 )=-0( K) /FOCTK 

0(K+1 )= ( Rh ( K )-E{ K)*E(K) J/FOC7R 

NGR = NGR ID- 2 

DC 130 K— 1 » N GR 

J = NGRIC-K 

P(J) = A( J+ 1 )*P( J + 1 )+B( J+ l ) 
1F(P(J),LT.PVAP) P( J)=PVAP 
CONT 1NUE 

w'R I TE ( 6 , 1 35 ) < P< J ) »J=1 » NG RI D ) 

FORM AT ( 2 X , 1 4E9 .2 ) 


INTEGRATE PRESSURE DISTRIBUTION TO GET X AND Y FORCE COMPONENTS 


UC 1 40 
ARG1 ( K ) 
AKG2 ( K ) 
A 1 
A2 
6 1 
3 2 
DO 
B2 
B 1 
C 1 
C 2 
DO 


K = 1 » NGR I D 

= P ( K ) *CTH( K ) 

= P< K ) +STH ( K) 

= ARG1 ( 1 ) +APG1 ( NGR IC ) 
= ARG2 ( 1 ) +ARG2( NGR I C ) 
= 0.0 


= 0.0 

150 K= 2 * NGRD * 2 
= 62+ARG2 ( K ) 

= Bl-MRGKK) 

= 0.0 
= 0.0 

160 K = 3 * N GR » 2 
C 1 =C 1 +ARG 1 ( K J 
C2=C2+ARG2( K ) 

DTFET=2 . * OX/BO 

F ACTR=FLC AT ( NF 1 LM J * D T HE T/3 . 
XO=F ACTR * ( Al+4. *61+2 . *C 1 ) 
YQ=F ACTR* ( A2 +4 
IF ( NSOLN .EC . 1 ) 

IF ( NSOLN .EC . 2 ) 

XC=XQ*F ACTR 
Y Q=YQ*F ACTR 
F X = FX + XQ 
F Y=FY+YQ 


*62 + 2 . *C2 ) 

F ACTR=-BD*BL/2 
FACT R=»BD*BL / 2 . 


CONTINUE 

GO TO ( 50 0 . 5 20 . 54 0 . 56 0 • 530 ) . KGUNT 

F 1 =F X 

F2=FY 

I F ( ( KOF < . E 0 . 0 > . AND . ( KGFC.EQ .0 ) ) RETURN 

IF(KCFK.EO.O) K OUN T = 4 

IF(KCFK.EC.O) GO TO 60 

GO TC 590 

AK yy=( FX-F 1 ) /DELS 


uu 
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LIST UTILITY 

AKYX=< FY-F2) /DELS 
GC TO 590 

540 AKXY=( FX-F1 ) /DELS 
AKYY= ( FY-F2)/0£LS 

C WRITE (6,6571 AKXX , A K X Y , AK Y Y 

657 FORMAT (5X, 'BEARING S T IFF . ( S CUEEZ ) : ' ,3 ( 5X * 013,6) ) 

AK11=-AKXX 
AK22=-AK Y Y 

A K 1 2=“ 0»5*(AKXY+AKYX) 

WRITE (6,657 ) AK1 1 , AK 12 , AK22 
WRITE (6,658) FI ,F2 

658 FORMAT (5X, 'FORCES ON THE R C 7 OR I • , // , 3 ( 2 X , 0 1 3 . 5 ) ) 

IF( KOFC.EC.O ) RETURN 
GO TO 590 

560 ACXX=( FX-F1 ) /DELST 
A C Y X = ( FY-F2)/OELST 
GO TO 590 

580 AC XY= ( F X- F 1 ) /DELST 
A C Y Y= ( FY-F2 ) /DELST 
G C TO € 0 C 
590 KOUNT=KOUNT+ 1 
GC TO 60 
600 CONTINUE 

A C 1 1 =« AC XX 
AC22=-ACYY 
RETURN 

300 WRITE< 6 , 700) 

RETURN 

700 FORMAT ( 1 HI// 5X * FOUR I ER—SER I EE 2-D OPTION NOT READY FOR USE'//) 


//?? 

// 

//?? 


322 RECORDS PRINTED 


END OF LIST UTILITY 


APPENDIX B 


Simple System "Driver" Fortran Listing 
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LIST UTILITY 


//NASA JOE 04 1 30 » 1 7 6 £ 3 A C AM S • • . M S GLE VEL= ( 2 • 0 ) 

/*J03P ARM TIME-9 
// EXEC FORT 
C 

C FOUR DOF RCT OR/BEARING/STATOR SYSTEM 

C 

DIMENSION XS{ 10023 * Y S ( 1002) * X8S( 1 002 ).YdS(1002).XREL{1002)t 
IYREL ( 1 002 ) *NPTS( 2) . INC( 2) 

D I MENS IGN *T(2002),Y7(2002) 

D I MENS IGN LINTYP{2) * INTEQ( 2) 

CALL FLCTS 

CALL PLOT ( 1 .0 , 1 . 5 t-3 ) 

I N = 5 
IC = 6 

5 CONTINUE 
KUNT=0 

SEAD(INi IOJBOi dL *3C»VISC»TH1 »7H2»P81 ,PB2 ,E ,CPM1 , PHI 1 , PVAP 
10 rCRMAT ( 5 E 1 5 • 7 ) 

RE AC ( IN', 15)N6RID*NS0LN, NPO RT ,NC YC , NDTP C , KOFK , KOFC . NF I LM 
1 5 FORMAT (1615) 

P 1 = 3 1 1 A 1 £92654 
WRITE! 10,20) 

20 FORMAT ( 1 HI///) / 

WR ITE ( 13 , 25)ED, BL, BC, VI SC, TH 1 ,TH2 ,PB l ,PB 2, E, CFH, PH II, PVAP 
25 FORMAT! 2X12E10. 3///) 

W R I T E ( IO*30)NGRID»NSCLN* NP ORT ,NCYC,NOTPC, KOFK ,KOFC .NFILV 
30 F ORM AT (2X101 10) 

READ (IN, 1 0 ) R MA SS * SM A S S t RFO RCE ,SKX, SKY 
READ(1N,10)X,Y,XDT,YC7,XB,YE.XBT,Y8T,WX,WY 
RE AD( 1 N , 15 ) NPR I NT , KPLCT , LI NTF ,NSK IP , KLUE 
R M AS S= RM AS S / 386 • 

SM ASS= SMASS/ 586 • 

Phil =P HI14PI/180. 

T AU 1 =6 0 • /CPM 1 

OT =T AU1/FLCAT (NCTPC ) 

NTS=NC YONDTPO 1 
OM 1=PI +CFM1/30, 

WRITE! 10 ,35 ) 

35 FORMAT ( 1 HI //5X»NT • , E X , • T I ME • , 1 1 X , » X' , 1 1 X , ‘ Y • , 10X , * XS' .lOX.'YB' , 
1 10X,*RX',10X,*RY' , 1 OX, ' SX • , 1 OX , • SY ' // ) 

DO 50 NT = 1 , NT S 
NTM1=NT- 1 
T=DT*FLO AT(NT“1 > 

IF (NT . GE .NPR lNT>WRIT£(lO,45)NTMI,T,X,Y,XB,Ye,XDT,YDT,X3T,YBT 
I F( KPLCT .EG. 0 ) GO TO 43 
M A A= (NT- 1 )/NSKIP 
M AB=NSK IP*VA A 
M A C= NT — 1 

IF< (NT.EC.l) .OR. (MAB.EQ.MAC ) ) GO TO 40 
GO TO 43 
40 KLNT=KUNT+1 
XS (KUNT ) =X 
Y S(KUNT ) = Y 
X B S ( KUNT ) =XB 
YBS( KUNT )=YB 
43 CONTINUE 

45 FORMAT (2X, I5.9E12.3 ) 
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LIST UTILITY 

• » 


30 


100 


ISO 


0M1T=0M1*T 
A R G=C M IT 4 PHI 1 
SARG=SIN< ARG ) 

CAR G= COS( ARG ) 

CALL S CUE E Z ( BC . BL , BC , V I SC . T H 1 , T H 2 • PS l * PB2 , X 
. XQT »YBT» AKI 1 « AK1 2 » AK22 , AC 1 l ,AC2 2 »F 1 *F2 ♦ NGR 
! KPFC, KFILMi PVAP) 

Rf X=RFORCEtCARG+Fl 
RF Y = RF OR CE* S ARQ + p £ 

SFX=-F 1-SKXOXB 
SFY=-F2-SKY* YB 
R X A= RF X / RM AS S 
r,YA=RFY/RMASS 
SXA=SFX/SVASS 
5 Y A = EFY/SMA5S 

xci=xot+rxa*ot 

YCT= YDHRYAAOT 
XE7=XDT+SXA*QT 


_ Y.XOT, YDT.XB.YB, 

NGR ID » NSOLN» NPORT iKOFK * 


i tw* 
+ *Y 


YGT= YB 
X= X + XD 


+S YA*DT 
*DT 


50 



)=LINTP 
'KUNT 


Y=Y+YDT ADT 
XB=XB+X9TADT 
YE=YB+YBT*0T 
CCNT INUE 

IFIKFLCT.EQ. 0) GO TO 100 
INTECtli 
INTEC. _ 

L I NT YP ( 

L INTYPI 

NPTS ( 1 

NPTS (2 )=KUNT 
INC( 1 ) =1 
I NC ( 2)=t 
DO 80 K 3 1 » KU NT 
XT(K)=XS(K) 

YT(K)=YS(K) 

XT (KUNT+K )=X9S( K ) 

YT(KUNT + K»=Y13S(K) 

XREL(K)- XS(K)«XBS(K) 

YREL < K J= YSIK)-YBS(K) 

CALL PLOT! 12 .0,0, 0,-2 ) 

CALL GRAPHIC AS, YS.XES. Y3S, XT ,YT .NPTS, l NC , L I NT YP , ( NTEQ , KUNT ) 
CALL PLOT( 12 ,0.0 ,0,-2 ) 

NPST=KLNT 
1 NK = 1 

I NTQ=69 

CALL GRAPH2 ( XREL , YRcL »NPST , INK, LINTP, INTO, KUNT ,0C j 

CONT INUE 

1 F ( KLUE . EG ,0 1 GO TC ISO 
GO TO 5 

CALL PLOT (10.0,0.0,999) 

STOP 

END 


SUBROUTINE GRAPH l(XS,YS,XBS.Y3S, XT, YT, NPTS , INC ,LINTYP, l NTEQ , KUNT ) 
DIMENSION XS{ 10 0 2) ,YS( 10 02 1 , >BS( 1 002) ,Y8S( 1 002) »XREL< 1002 ) , 
lYRELllOOe ) » N P T S ( 2 ) , l N C ( 2) 

DIMENSION LINTYP(2 ) , INTEQ( 2) 

DIMENSION XT ( 20C2 ) , YT ( 2002 | 
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L I ST UTILITY 


10 


20 


4 0 


7b2 


<2UNT=2*KUNT 
K21NTI = K 2 UNT + 1 
K2UNT2=K2UNT+2 
<LNT2=KUNT+2 
KLNT 1= KUN T+ 1 

CALL SCALE ( XT , 7 . 5, K2UNT , I NC ( 1 > ) 
CALL SCALE ( YT * 7 . 5 , K2UNT , I NC ( 1 ) ) 

X S < K LN T 1 )=XT< K2UNT1 ) 

YStKUKTl ) = Y T ( K2UNT 1 ) 

XS1KUNT2 ) = XT (K2UNT2 ) 

YS ( KUNT2 > = Y T ( K2UNT2 ) 

IF (XS( KUNT2) .GE .YS( KUNT2) ) CC TO 10 
X B S ( KUNT 2 )=YS< KUNT2) 

YES( KUNT2 )=YS( KUNT2 ) 
XS(KLNT2)=YS(KUNT2) 

GO TO 20 

Y S ( KUNT2 ) = XS ( KUN T 2 ) 

X8S( KUNT 2 ) = XS( KUNT2 ) 

YBS( KU NT £ )-XS{ KUNT 2 ) 

XEStKUNTl )-XS( KLNT 1 ) 

YE S ( KLNT 1 ) = Y S< KLNT 1 ) 


CONTINUE 
C ALL AXIS 


( 0 . 0 * 0 . 0 . • X-D ISPL * .“6 ,7,5,0. 0.XSCKUNT1 ) ,X5(KUNT2) > 
AXIS (0.0, 0.0, 'Y-DI SPL • , b , 7 . 5 , 90 . 0 , YS < KU NT 1 ) ,YSl KUNT 2) ) 
L I NE ( XS , Y S , NPT S ( 1) ,INC( 1) . LI NTYP( 1 ) , INTEQ ( 1 ) ) 

LINE ( XB£»Y3S»NFTS( 2) , INC (2) ,LINTYP( 2) . I NT EQ (2) ) 


C ALL 
CALL 
CALL 
CALL 

R E TURN 
c NC 

SUBROUTINE GRAPH2 < XR EL , YREL , NPST » 1 NK , L I N TP , I NT 0 , K UNT 
0 I ME NS ICN XREL( 1002 ) , YRELl 1002),XC(4 02),YC(402) 


SYMBOL ( 0 . 5, 7 • 7 E » 0 . T4 , ’ ROTOR AND STATOR OR 31 T 5 ’ , 0 . 0 . 23 ) 


, BC ) 


KUNT 1=KUNT+1 
KLNT2=KUNT+2 
DC 7 fc 2 K = 1 ,400 


Z = K/E0 . 

XC( K )=BC*COS( Z ) 
YC (K )=ec*s IN (Z ) 


C C NT INUE 

CALL SCALE ( XC , 7 . 5 , 4 0 C , 1 ) 

YC ( 401 ) = XC< 401 ) 

YC ( 402 ) = XC( 402) 

CALL AXIS(0. 0,0.0 
CALL A X I S ( 0.0, 0.0 
CALL L INE ( XC ,YC,400 , 1 
XREL ( KLNT 1 ) = XC ( 40 1) 

Y R E L ( K UNT 1 )=XC(401 ) 

XREL (KLNT2)=XC( 402) 

YREL(KUNT2 )=XC(4 02 ) 

C ALL L INE ( XREL, YREL , NPS T , I NK ,L_I NT,P, INTO ) 
CALL SYMBOL (0.5, 


•REL X-DISfL' ,-11 , 7.5, O.O.XCC401 ) ,XC( 402) > 
' REL Y—D ISPL' ,1 1 ,7.5,90#0,YC(401 ) »YC(402) ) 
0 , 0 ) 


5.7. 75,0. 14 ,* ROTOR ORBIT RELATIVE TO STATOR*,0 


i 20 ) 
RETURN 
END 


0 , 


D,L,C,y,Pj * * 9^ » 0 .j 

Bearing 


Damper 

One-Time Input 

El ement 


X R’ X S’ Y R ,Y S 

Xr,^s»Yr,X s 

Calling Mode from 
F.E. Program 

F X’ F v 


3F. 3F . 

[5^], [-4] 


J 3Xj 


Fig.l Input/Output of Damper Pilot Code 
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Fig. 2 Inertial coordinates 



Fig. 3 A typical aircraft engine damper configuration 


4 TO 8 IN. 



O-RING SEALS 



Fig. 4(a) Configuration frequently used in military applications 
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PISTON -RING 
SEALS 



Fig. 4(b) Configuration frequently used in commercial applications 
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Fig. 5 Unwrapped squeeze-film pressure solution domains for configurations 
shown in Figure 4. 


55 
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Figure 6 



LONG -BIDING '..CuUMCN 



(b) e/c=0.20 


Pressure distribution in circumferential direction and time 
of one cycle of circular orbit(long-bearing solution). 
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Figure 6 
( Cont ' d 



(c) e/c=0.60 



(d) e/c=0.95 


Pressure distribution in circumferential direction and time 
) of one cycle of circular orbit (long-bearing solution). 


pmak 
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(a) e/c=0.05 



Figure 7 Pressure distribution in circumferential direction and time 
of one cycle of circular orbit (short-bearing solution). 



mm fsi ^v> prmx-si 
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(c) e/c=0.60 


SMORT-atftRlNG SOLUTION 



(d) e/c=0.95 


Figure 7 Pressure distribution in circumferential direction and time 
(Cont'd) of one cycle of circular orbit (short-bearing solution). 
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Y 



Figure 8 Simple 2-mass , 4-degree of freedom. Test case 
(Same damper parameters as on page 23) 
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F i 9 • 9 



* ( AMs ) 

1 -■** & tg go 

(a) Rotor and stator orbits 



(b) Rotor orbit relative to stator 
(clearance circle shown) 


Nonlinear dynamic transient of simple 4 DOE sys tern ( See . F i g . 8 ) 
| F | = 1 00 lbs',«' = 150 rad/sec, M1 = M2 = 500 1 bs ,Kx=Ky=l 16000 lbs/in. 
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(a) Rotor and stator orbits 



(b) Rotor orbit relative to stator 
(clearance circle shown) 


Fig. 10 Nonlinear dynamic transient of simple 4 DOE sy s t em ( See . F i g . 8 ) 
| F | =200 lbs’,^ = 150 rad/sec ,M1 = M2=500 1 bs , Kx= Ky= 1 1 6000 Ibs/in. 
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(b) Rotor orbit relative to stator 
(clearance circle shown) 


Fig. 11 Nonlinear dynamic transient of simple 4 DOE system(See.Fig.8) 
| F | =300 1 bs',^= 1 50 rad/sec, M1 = M2=500 1 bs , Kx = Ky= 1 1 6000 Ibs/in. 
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(a) Rotor and stator orbits 



(b) Rotor orbit relative to stator 
(clearance circle shown) 


Fig. 12 Nonlinear dynamic transient of simple 4 DOE sys tem ( See . F i g . 8 ) 
| F | =500 lbs’,^ = 150 rad/sec ,M1 = M2=500 1 bs , Kx=Ky=l 16000 Ibs/in. 
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„ ( ) 

*b T ^ t j T — 6 

(a ) Rotor and stator orbits 



(b) Rotor orbit relative to stator 
(clearance circle shown) 


Fig. 13 Nonlinear dynamic transient of simple 4 DOE system(See.Fig.S) 
j F| =1 000 1 bs',^ = 150 rad/sec, M1 = M 2 = 500 1 bs , Kx= Ky= 1 1 6000 Ibs/in. 



FIG. 14 (OVERALL ARCHITECTURE OF BEARING IMPLENT CODE) 



Fig. 15 ARCHITECTURE OF CORE PROGRAM OF "BEARING ELEMENT" IMPLANT 
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FIG. 16 (OVERALL FE CODE ARCHITECTURE) 




FIG. 16 (CONT'D) 
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FIG. 16 (CONT'D) 
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2) 
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CONTACT UPDATES 

4) 

CHOICE OF ALGORITHMS 
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FIG. 16 (CONT'D) 



Fig. 17 (OVERALL PROGRAM FLOW) 
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